clear; clc;

b = 2*pi; c = 11.5;
x0 = 0.4; y0 = 0.8; z0 = 0.6;
delta = 1e-6;
N = 600;
transient = 50;

a_list = linspace(0.01, 5, 100);
LEs_all = zeros(length(a_list), 3);

parfor i = 1:length(a_list)
    a = a_list(i);
    LEs_all(i,:) = LEs(a, b, c, x0, y0, z0, N, transient, delta);
end

figure;
hold on;
plot(a_list, LEs_all(:,1), 'r-', 'LineWidth', 1.5);
plot(a_list, LEs_all(:,2), 'b-', 'LineWidth', 1.5);
plot(a_list, LEs_all(:,3), 'g-', 'LineWidth', 1.5);
xlabel('a', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('Lyapunov Exponents', 'Interpreter', 'latex', 'FontSize', 14);
legend('LE1', 'LE2', 'LE3');
title('Lyapunov Exponents of 3D-SIMM versus parameter a');
grid on;
hold off;

function LE = LEs(a, b, c, x0, y0, z0, N, transient, delta)
    eps_val = 1e-6;
    X = [x0; y0; z0];
    Q = eye(3);
    sum_LE = zeros(3,1);

    for i = 1:N
        % 避免除零问题，调整接近零的状态值
        for k = 1:3
            if abs(X(k)) < eps_val
                X(k) = sign(X(k)) * eps_val;
            end
        end

        % 系统迭代更新
        [x_next, y_next, z_next] = SIMM(X(1), X(2), X(3), a, b, c);
        X_next = [x_next; y_next; z_next];

        % 计算当前点的雅可比矩阵
        J = Jaco(X(1), X(2), X(3), a, b, c, delta);

        % QR 分解更新正交扰动矩阵
        Z = J * Q;
        [Q, R] = qr(Z);

        r_diag = diag(R);
        r_diag(abs(r_diag) < 1e-12) = 1e-12;  % 防止对数计算数值不稳定

        % 丢弃 transient 期，统计对数累积
        if i > transient
            sum_LE = sum_LE + log(abs(r_diag));
        end

        % 进入下一步状态
        X = X_next;
    end

    % 计算平均李雅普诺夫指数
    LE = sum_LE' / (N - transient);
end


function J = Jaco(x, y, z, a, b, c, delta)
    f0 = zeros(3,1);
    [f0(1), f0(2), f0(3)] = SIMM(x, y, z, a, b, c);
    J = zeros(3,3);
    for j = 1:3
        x_perturb = [x; y; z];
        x_perturb(j) = x_perturb(j) + delta;
        [fp1, fp2, fp3] = SIMM(x_perturb(1), x_perturb(2), x_perturb(3), a, b, c);
        J(:, j) = ([fp1; fp2; fp3] - f0) / delta;
    end
end
